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ABSTRACT 

We present the methods and preparatory work for our study of the collisional runaway 
scenario to form a very massive star (VMS, M* > 400 M 0 ) at the centre of a young, 
compact stellar cluster. In the first phase of the process, a very dense central core of 
massive stars (M* ~ 30 — 120 M 0 ) forms through mass segregation and gravothermal 
collapse. This leads to a collisional stage, likely to result in the formation of a VMS 
(itself a possible progenitor for an intermediate-mass black hole) through a runaway 
sequence of mergers between the massive stars. In this paper we present the runaway 
scenario in a general astrophysical context. We then explain the numerical method 
used to investigate it. Our approach is based on a Monte Carlo code to simulate the 
stellar dynamics of spherical star clusters using a very large number of particles (a few 
10 s to several 10 6 ). Finally, we report on test computations carried out to ensure that 
our implementation of the important physics is sound. In a second paper, we present 
results from more than 100 cluster simulations realized to determine the conditions 
leading to the collisional formation of a VMS and the characteristics of the runaway 
sequences. 

Key words: Galaxies: Nuclei — Galaxies: Starburst — Galaxies: Star Clusters — 
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1 INTRODUCTION 

Runaway collisions and mergers of massive stars following 
gravothermal contraction and core collapse in a young, dense 
star cluster provides a natural path to th e formation of a 
mass i ve object at the centre of the syste m febisuzaki et alJ 
l2QOlt [Portegies_Zwgxt_y^IhMjIlan||2nn^Ourkan. Freitag, & 
Rasio '21III ll hereafter GFR04). The basic idea goes back to 
the earliest studies of the quasar/AGN phenomeno n (e.g., 
ISoitzer fc SaslavJllQBfit IColgatel|l967t ISanderdll97fll . 

Runaways could easily occur in a variety of observed 
young star clusters such as the “young populous clusters” 
(like the Arches and Quintuplet clusters near our Galactic 
centre) and the “super star clusters” found in all starburst 
environments, including most galactic mergers llFiger et Id] 
Il999t iGallagher fc SmitJl99flh . The Pistol Star observed in 
the Quintuplet clust,ei~ ilFiger et alJll99Sh may well be the 
product of such a runaway, as demonstrated b y direct N- 
body simulations llPortegies Zwart fc McMillanll2002ll . If the 
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massive runaway collision product collapses to a black hole 
(BH), this scenario also provides a route to the formation of 
intermediate-mass black holes (IMBH) in star clusters. Dy¬ 
namical evidence for IMBHs at the centres of some globular 
clust e rs has been reported for many yea rs llGebhardt e t alJ 
l2002t iGerssen et alJl2002t Ivan der Marel el^ifll^Oodl . Ten- 
tative evidence has also been reported for an IMBH in the 
Galactic centre source IRS 13, whic h was recentl y r esolve d 
into a small cluster of bright stars llMaillard et al.ll2004fl . 
This could be the remnant of a much larger cluster that got 
tidally disrupted as its orbit aroun d the Galactic centre de¬ 
cayed through dynamical f riction llHansen fc Milosavlievid 
feoollGIirkan fc Rasiol200fitl . The very bright ultra-luminous 
X-ray source (ULX) associated with the young star cluster 
MGG 11 in the starburst galaxy M8 2 cou l d also h ave been 
prod uced through runaway collisions llPortegies Zwart et alJ 
l2004fl . A similar process may be responsible for the forma¬ 
tion of IMBHs in larger clusters, such as proto-globular clus¬ 
ters, as well as seed BHs in proto-galactic nuclei. These can 
later grow through a variety of proces ses including gas ac¬ 
cretion, stellar captures , and m ergers jMrmahv et, alJll99ll : 
iFreitag fc BeiR |2002b^jYu_^^rremajnT|200^ Hansen & 

Milosavlievic l20o7lwUthe'7rToebll200^TBrandfordll2004h . 
































2 M. Freitag, F. A. Rasio and H. Baumgardt 


The questions we address here are also central to our 
understanding of low-frequency gravitational-wave (GW) 
sources and the development of data analysis and detec¬ 
tion strategies for these sources. In addition, our work may 
be viewed as mainly relevant to galactic astronomy, the re¬ 
sults could also have applications in extragalactic astron¬ 
omy and cosmology. The direct injection into the centre of a 
galaxy of many IMBHs produced by collisional runaways in 
nearby young star clusters provides an important new chan¬ 
nel for building u p the mass of a central supermassi ve BH 
through mergers (IPortegies Zwart fc McMillan! 12002 '). It is 
possible th at this process is still ongoing i n our own Galac¬ 
tic centre lHan^ iL_fc Milosaylievi mM iKim et aljEool 
iGiirkan *&~^isioll2005ih In contrast, minor mergers of galax¬ 
ies are unlikely to produce BH mergers, as the smaller BH 
will rarely experience enough dynamical friction to spiral in 
all the wa^^Uh^enL^^Mh^noi^nassi^^alM^Volon- 
teri et al. 120031) . These ideas are also of critical importance 
for the design and planning of LISA, since the inspiral of an 
IMBH into a SMBH provides one of the best sources of low- 
frequency GWs for direct stu dy of strong field gravit y with 
a SMC^hase^nterferometerJCrAle^^^niornd|200^Phin- 
nev l200^Goihn^7^^iriiem20(^^^Siein2nn3^ MHiough 
the SMBHs found in bright quasars and many nearby galac¬ 
tic n uclei are thought to have grown mainly by gas accretion 
(e.g., [Sohajj|l - 982y H aehne lt et al ] 19981 ‘Fabian fc Iwasawal 

ll999tT^cheton^T20n3i ~ current models suggest that LISA 

will probe most efficiently a cosmological massive BH pop¬ 
ulatio n of lower mass, which is largely undetected llMenoul 
120031) . LISA will measure their masses with exquisite accu¬ 
racy, and their mass spectrum will constrain formation sce¬ 
narios for high-redshift, low-mass galaxies and, more gener¬ 
ally, hierarchicannodel^^fJm^AormatioiWejgj^Haehnelt 
& KauffmamJgOn^ ^^uriie^^Hoi^n^nn^Tvoionteri et alJ 
l2003t ISesana et al.ll2004i~ 

GFR04 concentrated on the early dynamical evolu¬ 
tion of young, dense star clusters. They performed dynam¬ 
ical Monte Carlo simulations for systems containing up to 
10 7 stars, and followed the rapid mass segregation of mas¬ 
sive main-sequence (MS) stars and the development of the 
Spitzer instability. They showed that, with a realistic initial 
mass function (IMF), these systems can evolve to core col¬ 
lapse in a small fraction of the initial half-mass relaxation 
time. If the core-collapse time is less than the lifetime of the 
most massive MS stars, all stars in the collapsing core may 
then undergo runaway collisions. The study in GFR04 was 
limited to the first step in this process, up to the occurrence 
of core collapse. About 100 simulations were performed for 
clusters with a wide variety of initial conditions, varying sys¬ 
tematically the cluster density profile, stellar IMF, and the 
number of stars. GFR04’s results confirmed that, for clus¬ 
ters with a moderate initial central concentration and any 
realistic IMF, the ratio of core-collapse time to initial half¬ 
mass relaxation time is typically ~ 0.1, in agreement with 
previous calculations. It was also found that, for all realis¬ 
tic initial conditions, the mass of the collapsing core (at the 
onset of collapse) is always close to ~ 10~ 3 of the total clus¬ 
ter mass, very similar to the observed correlation between 
central B H mass and total cluster mass in a variety of envi- 
ronments l|Ma gorrian et allll998t iMerritt fc Ferraresell200ll 

lHaring fc R.ixll2004lh 

In this and a following paper llFreitag. Giirkan. fc Rasiol 


l2005t hereafter Paper II), we go a step further in our study 
of runaways, by modelling the actual stellar collisions and 
following the early growth of the massive runaway product. 

This paper is organised as follows. In Section 2, we ex¬ 
plain in more detail the scenario for forming a massive object 
through runaway collisions and review the previous works on 
the subject. In Section 3, we present the numerical method 
and physical ingredients used in our simulations. Test sim¬ 
ulations are presented in Section 4. Finally, in Section 5, 
we summarise these results and introduce Paper II. In the 
latter, we will present the results of more than 100 simula¬ 
tions carried out to study runaway collisions in a variety of 
clusters. 


2 THE COLLISIONAL RUNAWAY ROUTE 
2.1 Important quantities 


Before reviewing the collisional runaway scenario, it is useful 
to define some quantities which are often referred to. 

If a star with mass Mi and radius Ri travels with rel¬ 
ative velocity across a field of stars of mass M 2 and 
radius R 2 and number density n 2 , the collision probability 
per unit time for this star, i.e., the reciprocal of its collision 
time, is 


^ 2 j — *Scoll V re l Tl 2 with 
^coh 


Scon = n(Ri + R 2) 2 


2G(Mi+M 2 ) \ 
(Ri + R2)(V%y ) ' 


(1) 


An important velocity scale for collisions between such stars 
is given by 


V 2 = 2G 


Mi + M2 

Ri -f- R 2 


= (617.5 kms -1 ) 


2 Mi + M 2 
M 0 


R@ 

Ri + R 2 


( 2 ) 


In most situations, IMj <C V, and the collision cross 
section is dominated by gravitational focusing, S'coii ~ 
2nG(Mi + M 2 )(Ri + i? 2 ). In a system where all stars have 
mass M*, radius R *, density n and a Maxwellian velocity 
distribution with 1-D velocity dispersion a v (<g V*), the col¬ 
lision time (after which, each sta r, on average, would have 
experienced one collision) is then llBinnev fc Tremainelll987l 
Eq. 8-125) 


f co ii ~ 2.1 x 10 12 yr 


10 6 pc 3 a v Rq Mq 
n 30kms -1 R* M* 


(3) 


Two-body relaxation plays a central role in the evolu¬ 
tion of most clusters considered here. The local relaxation 
time is llSpitzeilll987ll 


t r i x = 0.339 


G 2 n(M t y 


~ 4.8 x 10 yr x 


10 1 

( 1 

, 3 10 6 pc -3 , 


ln(7 c AL) 

V 30 km s 1 ) 

n 

V M 0 ) 


(4) 


where (M*) is the average stellar mass and N * the total num¬ 
ber of stars. For systems with a broad stellar mass spectrum, 
we use 7c = 0.01 in the Coulomb logarithm (see Sec. rom . 
The relaxation time defined this way only has a direct mean¬ 
ing for a single-mass population. In case of a mass spectrum, 
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it serves as a reference time but relaxational evolution may 
(and does) happen on a small fraction of f r ix- 

The core of a cluster is the central region where density 
and velocity dis persion are ap proximately constant. We use 
the definition of lSpitzeil Il987t Eq. 1-34) for the core radius, 

-Rcore = [9 ct„ 2 c / (47rGp c )] 1/2 (5) 

where p = (M*)n and underscore ’c’ indicates central values. 

2.2 Basic scenario: summary and expectations 

For the reader’s convenience we summarise here the colli- 
sional runaway scenario presented in GFR04 for the forma¬ 
tion of an intermediate-mass black hole (IMBH, 100 M@ < 
Mbh < 10 5 M@) at the centre of a dense stellar cluster. 
In the next subsection we briefly review previous works that 
have led to the formulation of this scenario or have pioneered 
its investigation. 

The basic idea is to create, through a sequence of colli¬ 
sions, a stellar object much more massive than what normal 
star formation can produce, with mass of a few hundreds 
to a few thousands Mg. We refer to these objects as “very 
massive star” (VMS). If such a star does not lose too much 
of its mass to winds by the time it reaches the end of its life, 
it is likely to undergo comp lete collapse into a blac k hole 
(BH)^withou^nas^jectioi^^T^^^Kaloger as Heger 
et al. l20n^^ ienc^pro^cin^u^^^^^^ 

Stellar collisions are extremely rare in most astrophysi- 
cal environments. They can only play a role in systems with 
a high stellar density such as self-gravitating dense clusters. 
When the most massive stars initially present in a cluster 
explode as supernovae, the cluster expands significantly as a 
result of the mass loss and collisions stop (if they were ever 
occurring). Therefore the formation of a VMS has to occur 
before massive stars evolve off the MS (including the giant 
phase would only make a small difference), i.e., for a clus¬ 
ter containing stars up to ~ 100 Mg, within a few million 
years. It is possible that a cluster could be bom in a colli- 
sional state, i.e., with an average collision time shorter than 
this, so that it would have been even more collisional while 
its stars were still on the pre-MS. This has been discussed 
as a way of creating all stars more massive than ~ 10 Mq 
iBonnell et al.lll998tlBallv fc ZinneckedlioOfil and references 
therein). However, we focus here on the simpler case of a 
gas-free cluster with all stars on the MS initially. 

Two-body relaxation provides a mechanism to increase 
the central density of a cluster and possibly drive it to a colli¬ 
sional state. Although this process of gravothermal core con¬ 
traction, eventually leading to core collapse, also operates in 
single-mass clusters, it is much more important for realistic 
clusters with a broad IMF. The most massive stars will then 
be subject to dynamical friction and drift to the centre. In 
all realistic cases, the cluster encounters the Spitzer insta¬ 
bility, meaning that this mass-segregation process leads to 
the formation of central core of massive stars that decouples 
as a self-gravitating system before the most massive stars 
can achieve energy equipartition with lighter objects. The 
system of massive stars then experiences core collapse on its 
own, very short, timescale. The high concentration of mas¬ 
sive stars at the centre of the cluster eventually leads to a 
high collision rate. 


The process of rapid mass segregation and core collapse 
in clusters with a broad IMF was the subject of GFR04. 
There, the relaxation-driven evolution of clusters with a va¬ 
riety of structures and IMF was followed. A key finding of 
this work is that for a given cluster structure but for all 
power-law mass spec tra, dN„/dM„ oc M ,~° with a £ [1.4, 31 
as we ll as for Kroupa llKroupa et al.ll993frl an d lMiller fc Scalol 
<197Slf IMFs, the core-collapse time t c c , expressed in units of 
the initial half-mass or central relaxation time (t r h (0), t rc (0), 
see GFR04) only depends on the ratio of the maximum to 
the average stellar mass, p = M*, max /(M*). Furthermore, 
for p > 50, a regime reached by all realistic IMF, the de¬ 
pendency flattens to a constant t cc /trh(0) ~ 0.07 — 0.08 for 
Plummer models. Even more interestingly, this value, when 
expressed in units of f rc (0) appears to be independent of the 
cluster’s structure 1 , t C c/t IC (0) ~ 0.15. GFR04 also found 
that, in contrast to core collapse in single-mass clusters, the 
mass of the contracting core, instead of decreasing to zero, 
reaches a finite value, representing in all cases a fraction 
1.5 — 3 x IGF 3 of the total cluster mass. 

These findings led us to the following two simple expec¬ 
tations: 

(i) If the dynamical role of binaries can be neglected 
(see discussion in Appendix E), any cluster where the core¬ 
collapse time, tec — 0.15 trc (0) (for realistically broad IMF) 
is shorter than the stellar evolution time of the most mas¬ 
sive stars present, t* ~ 3 Myr (if the mass function ex¬ 
tends to ~ 100 Mq) will enter a phase of rapid collisions 
between the massive stars that drive the core collapse. If 
the central velocity dispersion at that time is not in excess 
of ~ 1000km s -1 , collisions are not too disruptive and can 
lead to growth by mergers (lLai et alJll993t iFreitae fc Bend 
l2005lf . Since the most massive object will have the largest 
cross section for further collisions, it is expected to grow in 
a runaway fashion, i.e., much faster than any other star (see 
the simple mathematical model in Section EH1 ). 

(ii) When a runaway occurs, the final mass attained by 
the central VMS cannot exceed about 3 x 10 -3 of the to¬ 
tal cluster mass. This is also in agreement with the mass 
of the final collision product f ormed in several previous N- 
body simulation s of runaways iPortegie s Zwart et al.lll999al : 

IPortegies Zwart fc McMillanll2f)02l) . 

In paper II, we put these expectations to the test through 
cluster simulations that include stellar collisions. We shall 
see that, by and large, our results confirm point |(i)| For 
large cluster masses (> 10' Mq), cluster evolution may be 
driven, from the beginning, by collisions, thus accelerating 
collapse so that runaway may happen at lower densities than 
predicted by |(i)| W e do not find support in our simulations 
for prediction |(ii) | but likely because our simulation method 
may not be able to predict correctly the final mass attained 
at the end of the runaway. 

An obvious and legitimate question is whether the con¬ 
dition tec = t* = 3 Myr is ever fulfilled in real clusters. This 
is illustrated in Fig. 0 where we show this condition for a 
variety of King models in a plane representing the (initial) 
mass M c i(0) and half-mass radius i?h(0) of the cluster. We 

1 Provided one can defin e a non-zero t rc (0), which is not the c ase 
for 7 -models with 7 < 2 iDehnenlll993t iTrcmaine et al.lll994l) . 
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Figure 1. Conditions for rapid core collapse. This diagram shows which cluster masses, radii and initial concentrations will lead to core 
collapse in less than 3Myr, i.e., before the most massive stars evolve off the MS (a necessary condition for collisional runaway to occur). 
We consider clusters that have initially the structure of King models, with various concentrations, parametrised by the dimensionless 
central potential, Wo. The IMF is assumed to be Salpeter between 0.2 and 120Mq ((M*) ~ 0.69 Mq). In GFR04 we have showed 
that, for this or other similar IMF, core collapse occurs in 0.15i r c(0), independent of Wo. Each solid line corresponds to the condition 
t cc = 3Myr for a given value of Wo (see labels on the right side of the frame). Below this line, the core collapse time is shorter, above, it 
is longer. The dots on the lines correspond to models with 10 000 stars in their (initial) core. For this IMF, only a fraction ~4x 10 —4 
of the stars are more massive than 50 Mq, corresponding of 4 such stars in the core. The arrows show how much of a decrease in the 
total mass M c \ or half-mass radius R h leads to a shortening of t cc by a factor of 2. 

We also show the position in the (M c i,i?h) plane of a variety of observed clusters. The circles (in red in the on¬ 
line color version) are the Milky Way globular clusters from the compilation by Harris ( [199.6 , updated on-l ine at 
http://physun.physics.mcmaster.ca/Globular.htmll; the (magenta) pentagons are LMC clusters iMackev &; Gilmore!l2003) . Stars 
(in orange) represent p opulous you ng clusters, “super star clusters” and the cluster G1 of M31. Da ta for the Arches and Quintuplet 
clusters ar e taken from | F igei] (|20Q4t). for NGC 1705-1 and NGC 1569-A from lHo &; Filippenkol Jl996l) and for MGG-9 and MGG-11 (in 
M82) from lMcCradv etaL |j2003|) . The s haded ellipse (i n gree n) is the region in which the nuclei of dwarf ellipticals and bulgeless spiral 
galaxies are observed lGeh^^^al.ll2002j: IWalcher et alJl2005l) . For MW globular clusters, an age of 10 Gyr was assumed and the total 
mass was increased to correct for the decrease of the stellar average mass due to evolution of massive stars (but no account has been 
made of tidal stripping of stars). 
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have assumed a 0.2 — 120 Mg Salpeter IMF ((m) ~ 0.69 Mg, 
/x ~ 174). We chose M c i(0) and i?h(0) as parameters be¬ 
cause their present-day values are more easily accessible to 
observation than, say, the central density and because they 
probably vary less than other quantities during cluster evo¬ 
lution. On the other hand, these quantities do not by them¬ 
selves determine f rc ( 0 ) and, hence, t cc ; one needs to know 
how concentrated the cluster was initially (as measured, for 
instance, by the ratio of core radius to half-mass radius), 
a piece of information very poorly constrained by observa¬ 
tions of evolved clusters. For the sake of simplicity, we re¬ 
strict our discussion to the King family of cluster models, 
for which increasing Wo corresponds to increasing concen¬ 
tration. We note that, although most numerical simulations 
so far have been done with moderate initial concentration, 
Wo ^ 6 , observations of young clusters suggest they may 
initially have very small cores, correspo nding to Wo ^ 8 
ICampbeiTer al,lll992& [Moffat et alJl994 but.. se^McCrady 
et al. 12003 . concerning the core radius of R136). For initial 
Wo values ranging from 1 to 12 , we have plotted the line 
indicating clusters with t cc = 0.15t rc (0) = t* = 3Myr. Be¬ 
low the line, a cluster with that Wo would have shorter t cc - 
Points of various shapes denote observed clusters (young or 
old; see caption). We see that, even for a relatively moderate 
Wo = 7, a few clusters inhabit the region of parameter space 
for which one expects fast core collapse, possibly leading to 
collisional runaway. 

There is an important caveat to be made, concerning 
the interpretation of the results of GFR04 and Fig. 0 For a 
0.2— 120 Mg Salpeter IMF, only a fraction ~ 4x 10 ~ 4 of the 
stars are more massive than 50 Mg. To have 10 4 stars with ~ 
4 such massive stars within the core of a Wo = 3 or 8 cluster, 
the total number of stars must be 4.2 x 10 4 or 1.9 x 10 s , 
respectively. These values are indicated, for the various Wo, 
as black dots on the corresponding curves of Fig.Q Actually, 
our result that, for sufficiently large number of particles, the 
core collapse occurs on a given fraction (0.15) of the initial 
central relaxation time t rc ( 0 ) does not necessarily imply that 
only massive stars initially in the core (where the relaxation 
time is ~ t rc (0)) are responsible for the process. If it were 
the case, t cc would be of order of the dynamical friction 
timescale in the core, t d f iC — /i - 1 trc( 0 ) ^ 0 . 02 t rc ( 0 ) which is 
is much shorter than the value we find. This indicates that 
massive stars initially outside the core have time to reach 
the centre. 


2.3 Previous works on collisional runaway 

Colgate (1967) was the first to discuss the possibility of form¬ 
ing stars much more massive than initially present in a dense 
cluster through a sequence of stellar collisions. He suggested 
this mechanism as a way to create a large population of 
massive stars in a galactic nucleus to explain the quasar lu¬ 
minosities through enhanced super-nova rates. He pointed 
out that, provided the collisions always result in mergers 
with little mass loss, a runaway situation should ensue, with 
one star growing to a very large mass in a short time but 
also estimated that its growth would actually terminates at 
~ 50 Mg because, assuming a R, oc M* relation, the run¬ 
away object would then become too diffuse to stop 1 Mg 
impactors (“transparency” problem). A more realistic mass- 
radius relation for MS stars more massive than ~ 30 Mg (see 


Eq. 1121 below) actually corresponds to nearly constant pro¬ 
jected mass density M*/J ? 2 so that this argument does not 
apply if the growing star has time to contract back to normal 
MS structure between collisio ns. 

Following lColgatel lll967l) . one may gain some insight to 
the runaway mechanism by considering an idealised situa¬ 
tion in which one star of mass M(t) and radius R(t) grows 
by merging (without mass loss) with stars of mass and ra¬ 
dius m, r. One assumes M m (and R^$> r) and a constant 
density n of light stars, with a (3-D) velocity dispersion 0 - 3 . 
Then, using equation 0 , the growth rate of the massive 
star reads 


dM m , _\2 / 2G(m + M) 

= ~ mna 3 n (r + R 1 + , \ 2 ; 

dt tcoll V V + R ) a 3 


„ _ 1 7 , r r. Mo ( M 

c^2irGa 3 L nmMR = - ( - 

3 to \Mo 

with tg 1 = 2irGa 3 1 nmRo. 


1+/3 


( 6 ) 


This holds for strong gravitational focusing and a power-law 
mass-radius relation, R = Ro(MfMo ) 13 . The solution of this 
differential equation, for M(t = 0) = Mo is 

M(t) = - -—- -jjp with f div = (7) 

Hence, in this toy model, M becomes formally infinite after 
a finite time fdiv, if the exponent f3 is positive. More detailed 
analysis of the evolution of the whole distribution of stellar 
masses through use of the so-called “coagulation equation” 
also leads to the condition 0 > 0 for runaway to be possible 
jLeeil993l EoOOt Malvshkin fc GoodmanfeOOll). This kind of 
approach, ignoring stellar dynamical effects as it does, can 
only serve as a preliminary guide. A rough estimate of f d i v 
in a static cluster core would be 

= 2n(3GnmRo ~ (8) 

where t d yn is the core dynamical time. One sees that this 
is a very long timescale because R c /Ro is typically (much) 
kirgeiMhan 10 ®. 

Sanders ( 1970t) investigated the possibility of runaway 
collisions in dense galactic nuclei (without a central (I)MBH) 
using a more rehned model than Colgate’s. The evolution of 
a population of stars subject to collisions was followed using 
a “particles-in-a-box” Monte Carlo method. The outcome of 
the collisions (occurrence of merger, amount of mass and en¬ 
ergy loss) was deter mined using a generalisat ion of the semi- 
analytical method of lSpitzer fc Saslawl lS 196fih . The structure 
and dynamics of the cluster were not resolved. Instead, the 
system was treated as homogeneous within a spherical do¬ 
main, its size and density being evolved by considering the 
amount of mass and energy lost through evaporation of stars 
and collisions, respectively and assuming permanent virial 
equilibrium. All the gas lost in collisions was recycled into 
stars. The velocities were picked from a Maxwellian distri¬ 
bution with equipartition between stars of various masses. 
Thanks to a shallower M-R relation with /3 = 0.7, there was 
no transparency saturation to the growth of a VMS, which 
was followed from 0.5 Mg to more than 300 Mg in a 10 7 Mg 
nucleus 2 . A cluster model 10 times more massive, with the 


2 |Lightmari_^_Sha£ir^j fl978) have stated the transparency prob- 
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same initial velocity dispersion of ~ 500 km s -1 , did not ex¬ 
hibit any sign of runaway because the velocity dispersion 
raised above lOOOkms -1 and thus collisions became disrup- 

ht m a n & Shapiro (1978), citing unpublished work by 
Fall and Lightman, exposed the conditions for the onset of 
a collisional runaway using a simple evaporative analytical 
model for the contraction of the core of a single-mass globu¬ 
lar cluster. They found that the core must evolve to a state 
containing only a few hundreds to thousands of stars with a 
velocity dispersion of order 200 kms _1 , but such approach 
lacks physical ingredients such as mass segregation which is 
key in the evolution of more realistic systems. At any rate, 
an important point made in this work was that, in the col¬ 
lisional stage, the stars should experience mergers at such 
high a rate that they should have no time to recover ther¬ 
mal equilibrium between two collisions and could well stay 
bloated, hence br inging back the transparen cy problem. The 
famous paper of iBegelman fc Reed ( 11 978 fl introduced the 
process of runaway collisions to a larger audience as “one of 
the quickest routes to the formation of a massive object in 
a dense stellar system”. However, although they mentioned 
mass-segregation as a way to increase the density of massive 
stars on a shorter timescale, a self-consistent picture of the 
evolution a cluster subject to relaxation and collisions was 
still miss ing. 

Lee jl98 7) and lOuinlan fc Shapirol dl99(]l) were first to 
study the role of collisional mergers in numerical models self- 
consistently resolving the structure and evolution of stellar 
clusters (without a central massive BH). They applied very 
similar simulation methods and assumptions to quite dif¬ 
ferent systems. Using codes that solve the Fokker Planck 
(FP) equation, they had to keep the stellar mass function 
discretized into “components” and represented the cluster 
as a finite set of distribution functions (in energy-space), 
one for each stellar mass. This approach has the advantage 
of producing results virtually devoid of noise but imposes a 
very artificial treatment of stellar evolution and collisions. 
Stars in the same mass component have to share the ex¬ 
act same properties, including some average age updated as 
time passes and merger products are added to the compo- 


lcm in the following way. The growing star will be unable to stop 
a smaller impactor when the latter has more kinetic energy (at 
infinity) than what is required to punch a hole of its own cross 
section through the massive object, i.e., to unbound the mass it 
sweeps by passing through it. Neglecting the central mass concen¬ 
tration of the star, this leads to the (very approximate) condition 


M 2 GM 

lF r ~R 


ma\ , 


where M and R are the mass and radius of the runaway ob¬ 
ject, m and r the typical values for impactors, and <73 their 
velocity dispersion. Assuming a power-law mass-radius relation 
R* oc M? this translates into the following relation for the “sat¬ 
uration mass” of the runaway object, 

/»nw- 2 ) 

Afmax ~ —2 ' J 

with 7 ;) = Gm/r. A typical 0 value for MS stars (i.e., assuming 
the runaway object stays on the MS) is 0 ~ 0.5 and, except for 
very massive and compact clusters, 173 < 0.3u* ss 300kms -1 , so 

Afmax 100777. 


nents, a possible cause of rejuvenation because the authors 
assumed complete mixing of the stellar gas during collisions. 
Collisions were treated as mergers without any mass loss. 
The mass and orbital energy of mergers of any mass have 
to be cast into the predefined mass components. Another 
questionable aspect of FP simulations, when applied to col¬ 
lisional runaways, is the applicability of this formalism to 
mass components containing a very small number of stars 
(sometimes less than one). These models also included the 
dynamical formation of binaries through 3-body interactions 
and their subsequent hardening (and ejection) as a central 
source of energy capable of reversing core collapse and turn¬ 
ing off collisions in clusters with a relatively low number of 
stars. It has since been realised that there is a high probabil¬ 
ity for collisions to occur when these 3-body binaries interact 
with other stars, leading to a very significant red uction of 
the heatoi!^Uie^ 3 rovid^Chemnfn^^Miii|l 99 j^Fregeau 
et al. 120041 . 

Both lLeel dl987fl and lOuinlan fc Shapirol lll990l) started 
with clusters where all stars initially have th e same mas s 
(0.7 and 1M@, respectively). In the work of liLea dl987th 
who was studying globular cluster models, collisions are ac¬ 
tually tidal captures assumed to lead to quick merger. He 
concentrated on cases for which there was no real runaway 
growth or significant speed-up of core collapse due to mass- 
segregation. In all models but one, the core collapse was re¬ 
versed by heating due to 3-body binaries or by mass-loss due 
to stellar evolution of mergers. As stars more massive than 
22.4 Mq were not allowed, the simulation had to be stopped 
for t he only cas e in which conditions for the runaway were 
met. iLee 1 lll98 7) suggested that a “typical” proto-galactic 
nucleus, with 10 8 stars and a half-mass radius of ~ 0.4 pc 
would be subject to the “merger instability”. The goal of 
lOuinlan fc Shapirol dl990h was explicitly to look for the on¬ 
set of runaway collisions as a way to create a very massive 
star and, eventually, a seed IMBH that could grow into a 
massive black hole (MBH, Mbh > 10 5 Mg) at the ce ntre of 
a galactic nucleus, as suggested bv lBegelman fc Reesl dl97St) . 
The results of these simulations suggested that runaway col¬ 
lisions would occur provided that the half-mass relaxation 
time is shorter than ~ 10 8 yr (to beat stellar evolution) and 
IV* ^ 3x 10 6 (to avoid binary heating). The authors stressed 
that, as a result of mass segregation, the rise in the central 
velocity during collapse is only moderate and collisions do 
not become disruptive. Although not very realistic, these 
early studies made plausible the idea that successive colli¬ 
sions and mergers of main-sequence stars could lead to the 
formation of a ~ 10 2 — 10 3 Mg object. 

Through direct V-b ody simulations of clus t ers con tain¬ 
ing 2 000 to 65,000 stars. I Portegies Zwart et alJ lll999all and 
IPortegies Zwart fc McMUhmTfeOf^ showed that, in such 
low-TV* systems, dynamically formed binaries, far from pre¬ 
venting collisions (by heating the cluster and reversing col¬ 
lapse), actually enhance them by increasing the effective 
cross section. In these small systems, once the few massive 
stars have segregated to the centre, one of them will repeat¬ 
edly form a binary with another star and later collide with its 
companion when an interaction with a third star increases 
the binary’s eccentricity or induces a chaotic “resonant” in¬ 
teraction. The growth of this star is ultimately stopped by 
stellar evolution, or by the dissolution of the cluster in the 
tidal field of the parent galaxy. Given the small number of 
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stars in these simulations, the maximum mass of the colli¬ 
sion product is only ~ 200 Mq when mass loss from stellar 
winds is negligible. 

More recently, IPortegies Zwart et alJ f'2004) have im¬ 
proved on these early JV-body simulations, considering mod¬ 
els for tw o young clusters in the galaxy M82: MGG-9 and 
MGG-11 llMcCradv et al.l2003l) . Most calculations for MGG- 
11 were performed with 131 072 particles (“128k”), assuming 
a Salpeter IMF ranging from 1 to 100 Mq; those for MGG- 
9 used the same IMF and about 4 times more particles, 
in agreement with a higher estimated mass. For two simu¬ 
lations of MGG-11, the record-breaking number of 585 000 
particles was used in order to extend the IMF down to a 
more realistic 0.2 Mq (Salpeter) and 0.1 Mq (Kroupa). The 
initial conditions used were King models with dimension¬ 
less central potential ranging from Wo = 3 to Wo = 15. 
The authors found runaway growth of a VMS in all highly 
concentrated clusters (Wo ^ 9) with a half-mass dynamical 
friction timescale for 100 Mq stars shorter than 4 Myr. The 
mass reached was at least 800 Mq and up to 2700 Mq de¬ 
pending in the M-R relation. Incidentally they noted that 
no reasonable model for MGG-9 complies with these con¬ 
ditions but that MGG-11 may be in the right domain and 
suggested this may explain why an ultra-luminous X-ray 
source is (possibly) associated with the latter but not with 
the former. In contrast to what was found in smaller sys¬ 
tems, for models with > 10 5 particles (and without primor¬ 
dial binaries), a significant number of co llis ions (of order 
25-30 %) involve only single stars. IPortegies Zwart et alJ 
I2004I I also performed one (very computer-intensive) 128k 
simulation for MGG-11 with Wo = 12 and 10% primordial 
hard binaries and found a slight increase in the collision 
rate and final mass of the runaway object, suggesting that 
primordial binaries cannot prevent the runaway process by 
halting the core collapse process. Whether this holds in gen¬ 
eral, in particular for initially less concentrated clusters has 
yet to be established, through, e.g., a series of MC cluster 
simulations including primordial binaries (Fregeau et al., in 
preparation). 


3 SIMULATION METHODS AND PHYSICAL 
INGREDIENTS 

In the present work, we use a set of stellar dynamical simu¬ 
lations for collisional clusters to establish the conditions and 
manner under which a runaway occurs. Our basic numeri¬ 
cal tool is a Monte Carlo code similar to the one used for 
GFR04 but developed independently and presenting many 
differences in its structure. We describe this code briefly in 
Section cm The reason for using this different code here is 
that it was originally written to study high-density galac¬ 
tic nuclei including the effects of stellar collisions. In Sec¬ 
tion cm we explain our treatment of collisions and discuss 
various relevant aspects of their physics. 


3.1 The Monte Carlo code for cluster dynamics 

In the past few years, a new Monte Carlo (MC) code, 
ME(SSY)**2 (for “Monte Carlo Experiments with Spher¬ 
ically SYmmetric Stellar SYstems”) has been developed to 


follow the long term evol ution of dense clusters, with an em¬ 
phasis on galactic nuclei llFreitad200ltlFreitag fc Benzl200ll 
2002b). This code is based on the scheme first proposed by 
Henonl ll973ll ;o simulate globular clusters but, in addition 
to relaxation, it also includes collisions, stellar evolution, 
and, optionally, tidal disruptions and captures of stars by a 
central MBH through emission of gravitational waves. 

The MC technique assumes that the cluster is spheri¬ 
cally symmetric and represents it as a set of particles, each of 
which may be considered as a homogeneous spherical shell of 
stars sharing the same orbital and stellar properties. Unlike 
with the MC code used in GFR04, in the present imple¬ 
mentation, the number of particles may be lower than the 
number of stars in the simulated cluster (but the number of 
stars per particle has to be the same for each particle). An¬ 
other important assumption is that the system is always in 
dynamical equilibrium so that orbital timescales need not 
be resolved and the natural time step is a fraction of the 
relaxation (or collision) time. Instead of being determined 
by integration of its orbit, the position of a particle (i.e., 
the radius R of the shell) is picked up at random, with 
a R probability density that reflects the time spent at R: 
dP/dR oc 1 /Vi(R) where V r is the radial velocity. Unlike the 
code of GFR04, our scheme adopts time steps that are some 
small fraction / of the local relaxation (or collision) time: 
5t(R) ~ / min(T re i, Tcoii) with / of order or smaller than 
0.05 3 . Consequently the central parts of the cluster, where 
evolution is faster, are updated much more frequently than 
the outer parts. At each step, a pair of neighbouring parti¬ 
cles is selected randomly with probability P se i ec oc l/5t(R). 
This ensures that a particle stays for an average time St(R) 
at R before being updated. Because particles are evolved 
one pair at a time and to ensure perfect energy conserva¬ 
tion, the (spherical) potential produced by the collection of 
particles is represented by a binary tree structure allowing 
both determination of its value at any given R and its up¬ 
date after modification of the position or mass of a particle 
in 0(ln./Vp) operations, where N p is the particle number. 

The relaxation is treated as a diffusion process , with 
the classical Chandrasekh ar theory llGhandrasekharl 1196(1 
iBinnev fc Tremainelll987ll . similarly to what is done in the 
code used in GFR04. This treatment shares many assump¬ 
tions with the methods that are based on integrating the FP 
equation directly. Unlike those methods, ours is based on 
particles and hence allows us to incorporate further physics 
more naturally, e.g., a continuous mass spectrum and the 
inclusion of an anisotropic velocity distribution. The most 
important addition for the purposes of this paper is stellar 
collisions. Unlike diffusive relaxation, collisions cannot be 
treated as a continuous process but are individual events, 
each of which may importantly affect the orbit, the mass or 
the mere existence of a particle. When a pair is selected, one 
computes the collision probability between a star of the first 
particle and a star of the second, 

Ecoii = S co \\V r ^n5t. (9) 

Scoii is given by equation 0 and is connected to the max¬ 
imum impact parameter for collision, S'coii — ^max- Note 

3 One also makes sure that the time step is significantly shorter 
than any stellar evolution time (MS life-time) for stars around R. 
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that the total local stellar density n enters the relation, 
not m or ri 2 dFreitag fc BenJl2002bT ). A random number 
is picked from the interval [0,1 [ with uniform probability. If 
it is smaller than P co ii, a collision between the two particles 
has to be simulated 4 . The impact parameter is determined 
by b = & ma x\/A, where X is another random number from 
the interval [0,1 [. The other parameters specifying the col¬ 
lision, i.e., Mi, M 2 and are already known from the 
particles’ properties. The method for determining the out¬ 
come of the collisions is explained in the next section. 


3.2 Stellar collisions and stellar evolution 

In this section, we explain some important aspects of the 
“stellar micro-physics” of great importance for the runaway 
scenario in more detail: the outcome of collisions, stellar 
evolution and the interplay between them. We explain how 
we deal with these questions in our code and what are the 
associated uncertainties. 


3.2.1 General considerations about collisions 


Consider a collision between stars of masses Mi, M 2 and 
radii Ri, J? 2 - At large separation, their relative velocity is 
VX\ and impact parameter b. If they were point masses, 
their trajectories would be hyperbolas with separation and 
velocity at periastron, 


V m 


x + \/l + x 2 ’ 

KS (x + Vi + z 2 ) , 


( 10 ) 

( 11 ) 


with x = 



f Ri + R2 \ 

J 


If one neglects tidal effects until the stars touch, the relative 
velocity at contact is Kont = \J (KS) 2 + V* ■ Because of 
gravitational focusing, d m ; n is a more useful parameter than 
b to describe how central a collision is. 

For -C V*. gravitational focusing is important, 

Scon =! 2 tt(Pi + R2)G{M 1 1M 2 )(VXi)~ 2 and dP/d(d min ) = 


const, and most collisions result in merger with little mass 


loss, 8M/M <0.1 dBenz fc Hilldlll987 

1992llLai et al. 199:4 

Lombardi et alJ 2002: Sills et al.1 20021 

Freitae & Benzl200fill. 


When — few V*, a regime probably only reached in 
galactic nuclei in the vicinity of a MBH, the cross section 
is geometrical, Scon ~ n(Ri + R 2 ) 2 , most encounters have 
relatively large dmin (dP/d(d ul i n ) oc d m i n ) and are “fly-bys”, 
i.e., both stars survive and remain unbound. Only nearly 
head-on collisions are highly disruptive and may lead to de¬ 
struction of the smaller star or both l| Benz fc Hillslll987l 
Il992t lLai et~aT1ll99Bt iFreitag fc Bendl2002al l200Blb 

At low relative velocities, two stars can become bound 
to each other, i.e., form a binary, through dissipation of 


4 If each particle represents N s t stars, N s t identical collisions are 
assumed to take place. In this way, one may apply the result 
of the collision (new velocities, stellar masses,...) to the particles 
themselves. If the collision results in a merger, only one particle is 
retained in the simulation. If both stars are completely disrupted, 
both particles are removed. 


orbital e nergy in tides ne ar periastron, even if d m in > 
Ri + R 2 iFabian et m. For velocities typical of globu¬ 
lar clusters (10 — SOkms -1 ), such tid a l binaries form up to 
dmin /(Ri + RA ^ 2 — 3 dF’ortegies Zwart fc Meineiilll993t 
iKim fc Leelllifflffi b Their evolution is still a subject of de¬ 
bate and may be very complex, but given the fact that they 
form with very small pericentre separation and that the stars 
should sw ell due to conversion of tidally excited oscillations 
into heatjMd'^km ct ah i 19871 iKurnar fc Goodmanlll99 fii: 

lPodsiadlowskillll99fih . it seems likely that they will promptly 

merge. Consequently, we could try to account for tidal cap¬ 
tures by inc reasing th e effective stellar radii for collisions 
at low VXi keel 11987 1. but for simplicity and to stay on 
the conservative side when testing the runaway scenario, we 
decided to neglect this effect and only account for genuine 
collisions with d m i n < i?i + R 2 . 


3.2.2 Mass-radius relation 


One may expect the M-R relation to play an important 
role in determining collision rates through its influence on 
the cross section. Fig. 0 show various M-R relations from 
the literature. For stars between ~ 0.8 and ~ 150 Mq, we 
plot the radius at the beginning, middle and en d of the 
MS ([Schall^^tayl^g^jChaAoime^^^l^TOS^Lejeune 
& Schaerer T20tBli ^hiBi^^resenr^oA^Br^^ehaviour of 
the M-R relation beyond M* = 100 Mq is of particular 
concern because, within our assumptions, this will set the 
size of the runaway collision products and may have bearing 
on growth and rejuvenation rates. Hence, this may deter¬ 
mine when stellar evolution will terminate the growth and 
what mass the runaway star has attained at this stage. Un¬ 
fortunatel y, the radiu s of a VMS (M, > 100 Mq) is highly 
uncertain. Ilshii et "HI f l999l ~) have shown that a VMS may 
dev elop a very extended envelope of very low density (see 
also Stothers fc Chinlll997t iFiger et a I . 199? I IBaraffe et al.l 


12001 h Their 1000 Mq, Z = Zq model, for instance, has 
R t = 3000 Rq on the zero-age main sequence (ZAMS) but 
its envelope represent only 3% of the stellar mass so it is 
reasonable to assume this extended atmosphere will play a 
negligible role in collisions. Furthermore, this structure is 
due to opacity from metals and is not p resent if the metal- 
licity is low enough (IBaraffe et alJl2001i '). 

For the cluster simulations presented in this paper, 
we adopt a unique MS mass-radius relation constructed 
from Z = 10 -3 models from the following sources. Fo r 
M*/Mq < 0.4, we use Fig. 3 of lChabrier fc Baraffel ll200Cli . 
For the ranges 0. 4 ^ M t / Mq < 1 a nd 1 ^ M* /Mq < 120 we 
use m odels from lCharbonnel et alJ dl999ll and lSchaller et all 
lll992l) . respectively. Finally, for stars more massive tha n 
120 Mq, we apply a relation given by iBond et al.1 dl984ll . 


B - = L6S ®(^)“ 7 ' (12) 
Even though this expression is for stars more massive than 
10 4 Mq , it appears to match nicely with the M-R relation 
at M*/Mq < 120, as shown in Fig. 0 


5 Data available at http://obswww.unige.ch/~mowlavi/evol/ 
stev_database.html 
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Figure 2. M-R relations for MS stars 
from various authors. The M-R for low 
masses (0.01-1 Mq, in cyan in the colour 
version^i^fron^^ig^^^o^Chabrier Sz 
Baraffe <200(1 Z = Z G =0.02, age of 
5 x 10 9 years) For intermediate masses, 
we plot radii from the Geneva stellar evo¬ 
lution group for three metallicities: Z = 
0.0004,0.001,0.02 (in blue, magenta and 
red i n the onl i ne color version of this fig¬ 
ure ) iSchuller_et_a]Jl^9^ |Cha£bonn^ et al.l 
Il999t iLeieun^^^Schaereil 1200ll) . Short- 

dashed lines correspond to the ZAMS, 
solid lines to the radius when the star has 
lived half of its MS life-time (or 5 x 10 9 
years for low-mass stars) and long-dashed 
lines to the end of the MS phase or at 
an age of 10 10 years at low masses. The 
maximum mass considered in these series 
of models is 100 to 150 Mq. All other 
M—R r elations plotted here are for the 
ZAM S iBond^etaj Jl98^Jgtothers_jj Chinl 

lgQTyisM^t^al ToognBara^^^alll^OOir 

gchaereil|2Qflj). Star with higher metallic- 
ity have larger radii. At Zq, stars more 
massive than 100 Mq have a huge diffuse 
envelope, d ue to metal opacity. The re¬ 
lation from iBond Rt all ill 984 (R, ~ 
1.6 figfM./MQ) 0 ' 47 ) neglects this effect. 
It is established for M* ^ 10 4 Mq but 
matches nearly perfectly the models for 
M, ^ 100 M 0 . 


Focusing on low metallicity is reasonable: only if the 
metallicity is sufficiently lo w, is it c lear that a VMS will 
be stable while on the MS (lBa.ra.ffe et ~ai1l200lft 6 . that it 
will experience little evolution ary mass loss, and that i t will 
collapse into a BH as a whole dFrver fc Kalogerall2()0T ). 

We do not account for the increase of the stellar radius 
during the MS evolution and use the size of the star at half 
its MS lifetime, or at an MS age of 5 Gyr for stars with MS 
lifetime exceeding 10 Gyr. Although it would be more con¬ 
sistent to let the stellar radii evolve, the extra complication 
involved in such an improvement is not justified in face of the 
larger uncertainties introduced by other necessary simplifi¬ 
cations, in particular concerning the size and evolution of 
collision products (see below). By the nature of the scenario 
studied here, we follow only the first few million years of the 
cluster evolution. Hence, low-mass stars should be given a 
size smaller than their “half-MS” radius and closer to the 
ZAMS value. One can see in Fig. |3 that this only leads to 
a slight overestimate of the collision cross sections which, 
being dominated by gravitational focusing, are proportional 
to the stellar radii. More important is the case of the most 
massive stars, M* > 50 Mq, say, which are expected to dom¬ 
inate the collisional process and may significantly evolve over 
the period of time considered here. Their size increases sig- 


6 But the question of the structure, stability and evolution of a 
VMS subject to a steady bombardment of smaller stars is another, 
unsettled question. Collapse to an IMBH may occur before the 
runaway collision product has evolved to thermal equilibrium, in 
which case the stability issues may be less relevant. 


nfficantly during the late MS. Neglecting this may lead to 
an underestimate of their collision rates. A higher, possi¬ 
bly more realistic, size contrast between light and massive 
stars is likely to facilitate the runaway mechanism. However, 
by experimenting with large changes in the high-mass M~R 
relation, ranging from f?* = const to R * oc M*, we have 
checked that this is of relatively little importance. 

A more questionable simplification is to consider that 
all stars are initially on the MS. The runaway process has 
to start before the most massive stars in the IMF (M» ~ 
100 Mq) turn into compact remnants. This gives us at most 
3 Myr (see Fig. OK)- This is to be compared with the dura¬ 
tion of the pre-MS phase. After it has stopped accreting (and 
hence reached the “birth-line”), a pre-MS star less massive 
than ~ 6 — 8 Mq is several times larger than on the ZAMS. 
The time required for contraction onto the ZAMS strongly 
increases with decreasing mass; it is short er than 3 M yr only 
for stars more massive than 2 — 3Mq f Pallal 12002ft . Con¬ 
sequently, in a young cluster, assuming for simplicity that 
accretion ceased for all stars at the same instant, low-mass 
stars are more extended and less dense than high-mass ob¬ 
jects. During this phase their collision rate is higher than in 
the MS phase; however their encounters with more massive 
but more compact stars may result in the tidal disruption of 
these low-massol^ec^mmHiei^iai^WbleaiFMnergei^Zin- 
necker & Bate 12002ft . We leave these complications out of 
our present study. We consider the models presented here as 
the simplest possible ones that incorporate all the key phys¬ 
ical ingredients needed to investigate the runaway scenario. 
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Figure 3. (a) Main sequence lifetime for stars of various masses, 
(b) Relation between the mass of a MS star and the amount of 
helium produced in the core by hydrogen burning (at the termi¬ 
nal MS), expressed as a fraction of the total mass. Values for two 
metallicities are plotted. Data kindly provided by K. Belczynski. 
For masses larger than 120 Mq, extrapolation (linear in the quan¬ 
tities plotted) is used; the fractional amount of He produced is 
not allowed to exceed 0.56 iBond et al Jl984ll . 

More sophistication could be included in further works, if 
deemed necessary. 

3.2.3 Collisional rejuvenation 

The stellar evolution of collision products has only been ex¬ 
plored for t he case of low-ve l ocity mergers, relevant to globu¬ 
lar clusters IlSills et al ] 19971 .120011) . Such encounters are only 
mildly supersonic and entropy is nearly conserved. Hence, 
the structure of the merger can be established by sorting 
the mass elements from the parent s tars according to their 
entropy llLombardi et a aaa The main uncertainty 
about the evolution of these objects is the mechanism, if 
any, responsible for decreasing the amount of angular mo¬ 
mentum as the star relaxes to thermal equilibrium after the 
collision. For high velocity collisions, significant dissipation 
occurs and entropy sorting is questionable, not to mention 
that, in most cases, both stars survive the collision unbound 
to each other (“fly-bys”). 

In view of these difficulties, we used a very simple pro¬ 
cedure to set the stellar evolution of mergers, called minimal 
rejuvenation. We assume that, during a coalescence, the he¬ 
lium cores of both parent stars merge together, while the 


hydrogen envelopes combine to form the new envelope; no 
hydrogen is brought to the core. Furthermore, to assign an 
effective age to the merger, we assume that the mass of the 
helium core grows linearly with time during the MS and 
resort to stellar evolution models to provide the relation 
between the stellar mass and the helium core mass at the 
terminal age MS. This formalism is also applied to fly-bys 
during which part of the stellar envelopes is removed. In 
Fig- mi we plot the MS lifetime and the total mass of helium 
produced during the MS phase as func tions of the stellar 
mass (data provided by K. Belczynski; see lHurlev et all200(t 
IBelczvnski et alJ 1200211 . For simplicity and because the de¬ 
pendence on metallicity is weak, we use the Z = 10 -4 data 
for all simulations. In any case, the thermal timescale is al¬ 
ways assumed to be shorter than the average time between 
collisions so that the MS mass-radius relation is applied to 
collisions products. The validity of this hypothesis during 
the runaway growth a massive star through repeated merg¬ 
ers is questionable and is discussed in Paper II (Section 2.2 
in particular). 

3.2.4 Implementation and role of stellar evolution and 
mass loss 

ME(SSY)**2 implements a very simp le prescription for stel¬ 
lar evolution llFreitag fc Benzll2002bl) . While a star is on the 
MS, its mass and radius are kept constant. The duration 
of thf^MSj^^^j^riwiMa^^endlet^telkumTiodel^Hurley 
et al. 1200 0;'). It is plotted for two metallicities in panel (a) of 
Fig.mi For all simulations presented here, we use the data for 
metallicity Z = 10 -4 but it is obvious that the metallicity 
has little impact on Tms- On the other hand, the amount 
of mas^Joss^i^h^^dSMncreases^trongl^^itl^^^Vink 
et al. 1200 ill . While the possibility of runaway collisions in 
clusters with high metallicities, such as young populous or 
“super” clusters observed in our and other galaxies is cer¬ 
tainly of great interest, it is unlikely that a high -Z VMS 
will form an IMBH if it is left to evolve on the MS precisely 
because it should experience such high mass loss. Indeed, 
standard prescriptions for mass-loss rates indicate that stars 
more massive than ~ 120 Mq shed most of their mass on 
the MS if they are of solar metallicity. At Z ~ 4 x 10 —4 , 
~ 120 Mq stars lose only ~ 10 % of their mass but objects 
above ~ 500 Mq should evaporat e themselves nea rly com¬ 
pletely llLeieune fc Schaererll200il lKudritzkill2005) . There¬ 
fore, in this work, we only consider (very) low metallicity. 
Also recall that, in GFR04, we studied whether mass loss 
on the MS could decrease the binding energy of the cluster 
enough to reverse core collapse and found that, even for so¬ 
lar metallicity, this only happens in a very small domain of 
the parameter space, namely for clusters with a core collapse 
time already very close to the critical value of 3 Myr. 

We assume that all mass lost by the stars is expelled 
from the cluster. This is probably a reasonable simplifica¬ 
tion for clusters with a relatively low escape velocity, like 
globular clusters but a poor one for (proto-)galactic nuclei 
in which a significant fraction of the gas could be retained. 
This approximation is discussed in Paper II. 

In summary, for this work, the role of stellar evolution 
is only to set a clock against which core-collapse and colli¬ 
sions have to race. When core collapse requires more time 
than the MS lifetime of massive stars, it will be terminated 
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by SN mass loss (see Fig. 13 of GFR04). This also holds for 
clusters of very low metallicities because all stars with mass 
between ~ 8 M@ and ~ 40 M@ at the end of the nuclear 
burning p hases should explode as SNe, with considerable 
mass loss dFrver fc Kalogeral 12001 ). Stellar evolution of a 
VMS interrupts the runaway collision sequence. If it turns 
into an IMBH, it may continue to grow, although presum¬ 
ably at a smaller rate, by tidally disrupting the MS stars 
while they still dominate the central density. When all nor¬ 
mal MS stars more massive than about 20 — 25 Mg have 
turned into stellar BH with M > 10Mg, these objects will 
dominate the central region, expelling MS stars and thus 
quenching tidal disr uptions. iV-body simulation s of systems 
with TV* ^ 180 000 feaumgardt et al Jl2004albt) have shown 
that the central IMBH will likely form a binary with a com¬ 
pact remnant which, although not compact enough to merge 
over a Hubble time through emission of gravitational radia¬ 
tion, is probably very efficient at scattering off other stars, 
thus driving cluster re-expansion and preventing any signif¬ 
icant growth of the IMBH. These conclusions presumably 
do not apply to systems containing 10 6 stars or more in 
which the bina ry shoul d have a se paration small enough for 
quick merging feaumgardt et ahl predict a oc Re¬ 

cent simulations performed with a Monte Carlo code which 
incorporates primordial binaries have revealed the exciting 
possibility of the formation of two IMBHs. In models with 
more than 10 % primordial binaries, a first VMS starts grow¬ 
ing at some distance from the centre through binary inter¬ 
actions while the second forms in the centre, at the moment 
of core colla pse throug h c ollision between single stars, as 
studied here llGiirkan et al.lt2005 ). 

The kind of MC code we use in this work does not 
include binary dynamics and, more generally, is not suited 
for following the cluster evolution when driven by a very 
small number of central objects. Hence, we do not attempt 
to study the evolution of the cluster and its central object 
once one VMS has formed and reached the end of its MS 
phase. This would be better studied using IV-body methods 
to treat in detail the full dynamics of the central region. 


3.2.5 SPH collision simulations 

In the MC method, collisions are handled on an event-by¬ 
event basis, rather than through average rates and outcomes 
as in Fokker-Planck codes. To treat collisions with as much 
realism as possible, we can use the results of 3-D hydro sim¬ 
ulations performed wit h the Smoothed Particle Hy drody- 
namics (SPH) metho^jBen ^^QOyMonaghapJ^KiL ^asio 
& Lombardi ll99fitl . A set of some 14 000 SPH simulations 
of collisionsJoeBreer^MS^R^^va^oerfomiet^a^Freitag 
& Benz T^OOsfl ^^iu^ririna^luar was to include the ef¬ 
fects o f collisions in stellar d ynamical models of galactic 
nuclei ( iFreitag fc Benzll2002bl) . This focus determined the 
parameter range covered in this study: stellar masses from 
Afmin = 0.1 Mg to M max = 74.3 Mg, relative velocities in 
the range /V* ~ 0.03-30 and impact parameter cor¬ 
responding to d min /(Ri + R 2 ) = 0-0.9. Because it proved 
intractable to summarise the outcome of these simulations 
(stellar masses, orbital energy and deflection angle) through 
a set of fitting formulae, we implemented a scheme to inter¬ 
polate from SPH results in the four-dimensional parameter 
space (Mi, M2, Krf; d m in)- While this method proves ade¬ 


quate for galactic nuclei simulation in which velocities are 
high and it is therefore highly unlikely that mergers will 
leac^cMbnnatioi^^^tarmoi^jnassiv^^imWHjjj^^Fre- 
itag l200flFl^reito^^^l20^r it is not suited to the problem 
of runaway growth. 

I11 most cluster simulations reported here and in Pa¬ 
per II, we simply assume that all collisions result in merger 
with no mass-loss. We will see that this is a satisfying ap¬ 
proximation because, in all but the most extreme cases, the 
central velocity dispersion is and remains relatively small. 
To go one step beyond this zeroth order “sticky sphere” ap¬ 
proximation and test its validity, for a subset of runs, we 
use the SPH results to allow for non-merging collisions and 
collisional mass loss in the way described below. For graz¬ 
ing collisions not resulting in a merger, we do not account 
for the reduction of the relative velocity or non-Keplerian 
deflection angle. 

Complete disruption of both stars is extremely unlikely. 
Not only does it require a relative velocity of a few V* but 
also a nearly head-on geometry. Hence, in this work, we con¬ 
sider only two possible outcomes: merger and “fly-by”. Most 
low-velocity (I/j.“ < 0.1 V*, say) encounters result in mergers 
but, in our SPH work, we have not studied those in much 
detail. Such collisions require a large amount of computation 
time because, unless the collision is nearly head-on, the pair 
does not merge immediately after first contact but forms a 
bound binary which goes through many successive perias- 
tron passages until final coalescence. In most cases, the hy- 
drodynamical simulation was stopped before a merged star 
had formed. But it is very likely that any binary formed 
through a contact collision will eventually merge because 
the stars swell as their envelopes are shocked so that each 
pericentre passage is more dissipative. Hence, we can use 
our SPH results to establish the condition for merger. We 
have found the following parameterisation to correctly pre¬ 
dict the maximum merger impact parameter for all but a 
few SPH simulations: 


Amerg — Co ~b c\q + C 2 M 2 + c%l v -\- c^lvQ ~\~ c§l v M 2 ~\- 


with A = 


cqIZ + ciltq + cslZM 2 

drnii 


offi) , o(h) ’ q - M 1 /M 2 < 1 

ilj il2 


(13) 


V° 


and l v = log 10 -^L, V, (h) = 

V* 


2G(Mi + M 2 ) 

p (h) 1 p(h) 

| 1X2 


Index 1 indicates the less massive star; R^l are the radii 
enclosing half the mass of each star. The numerical coeffi¬ 
cient are (co, ci,... cs) = (0.525, —0.107, 6.84x 10 -3 , —2.03, 
0.525, 6.16 x 10“ 4 , 0.132, 0.526, -4.89 x 10“ 3 ). 

To determine the final stellar masses, we make use of 
the SPH mass loss results in the following way. First when 
M2 > M max = 74.3Mg, we try to rescale Mi and M2 
by some factor r/ bringing r/M 2 to 0.95 M max while keep¬ 
ing tjMi > M m in ■ This would conserve the mass ratio q but 
is not possible when q < q m i n = M m i n /M max ~ 10~ 3 . Such 
cases are relatively unimportant because, in typical cases, 
the VMS grows mostly by merging with 50 — 120 Mg stars 
so q > 0.025 even for My ms = 2000 Mg. We deal with 
them by rescaling Mi to 1.05M m i n , independently of M2. 
From the SPH results, we interpolate the fractional mass loss 
5 m = (<5Mi + <5M 2 )/(Mi +M 2 ) for the (Mi, M 2 , KS, dmin) 
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Table 1. Important quantities for the Plummer model. All quan¬ 
tities are in A-body units except times which are in FP units. 


Quantity 

Symbol 

Value 

Plummer scale 

R P 

3 7r/16 ~ 0.589 

Core radius 

Rcore 

0.417 

Half-mass radius 

R h 

0.769 

Central density 

Pc 

1.167 

Central ID velocity dispersion 

&v ,c 

0.532 

Mass within R c ore 

Mcore 

0.193 

Central relaxation time 

trc 

0.0437 

Half-mass relaxation time 

Ah 

0.0930 


parameters of t he collision (with Mi, M 2 possibly rescaled), 
as explained in iFreitag fc Bend (' 2005 ). Inspection of SPH 
results for fly-bys shows that only in very exceptional cases 
does the mass of any of the two stars increase. In fact, to a 
good approximation, we find the mass ratio to be conserved. 
We thus use this assumption to determine the individual 
masses from 8m- 


4 TEST SIMULATIONS 

Before embarking on the large-scale numerical exploration 
of the runaway scenario to which Paper II is devoted, we 
first demonstrate that our numerical tools are up to this 
task. Here we check that we can reliably model the following 
two key aspects of the scenario: (1) Fast segregation-driven 
core collapse in clusters with a broad mass function and 
(2) Influence of collisions in cluster evolution and onset of 
collisional runaway. 

4.1 Important quantities and units 

Keeping with the tradition, when no t stated other wise, we 
are using the A-body unit system llHeno defined 

by G = 1, M c i(0) = 1 (initial total cluster mass) and 
U c i(0) = —1/2 (initial cluster potential energy). As time 
unit, we prefer the “Fokker-Planck” time Tfp to the A-body 
unit Tnb because the former is a relaxation time while the 
latter is a dynamical time; they are related to each other by 
Tfp = A*/ ln( 7 c A*) Tnb. 

All test com putations performed here are based on the 
Plummer model IPlummerlllOllI : iBinnev fc Tremaine!11 987t) 
for which we recall some important quantities in Tabled 
We refer to GFR04 (Section 3 and Table 1) for more de¬ 
tailed explanations about units and the important physical 
parameters of a variety of clusters models. 

4.2 Core collapse without collisions 

4-2.1 Comparison with direct A -body simulations 

In IFreitag fc Bend (1200 ill , the ability of ME(SSY)**2 to 
follow the relaxational evolution of clusters has been suc¬ 
cessfully tested against other simulation methods for clus¬ 
ters with single-mass, 2-component or continuous mass func¬ 
tions. However, in that paper, no comparison was done 
with direct A-body results and the only continuous mass 
spectrum considered was a relatively narrow 0.1 — 1.5 Mq 



Figure 4. Core collapse of a single-mass Plummer model. We 
compare the results of ME(SSY)**2 with A p = 300 000 particles 
(solid lines) to those of a direct A-body run with A p = 65 536 
\ dotted lines, in magenta in the on-line colour version; from 
iBaumgardt et H.l.l toO.'li . We plot the evolution of various La¬ 
grange radii, i.e., radii of spheres enclosing the indicated fraction 
of the cluster mass. For the A-body simulation, only fractions 
larger or equal to 1 % are tracked. Fractions of 3, 4,. ..9% are 
tracked in the A-body run but not in the MC computation. A 
value of 7 C = 0.09 was used in the Coulomb logarithm when con¬ 
verting A-body time units into FP units to get the best agreement 
in core collapse times. To improve clarity, the curves have been 
smoothed using a sliding average with Gaussian kernel. 


Salpeter IMF (p = M*, max /(M») ~ 2.2). Given the im¬ 
portance of core collapse in clusters with p > 100 for the 
runaway scenario and the dearth of published data about 
this situation, we decided to carry out new comparisons be¬ 
tween ME(SSY)**2 and direct A-body integrations of the 
core-collapse evolution of clusters with various IMFs. 

The A-body simulations reported here were done by 
H. B. with the collisional Aarseth A-body code NBODY4 
ArmeUi^^) on th e GRAPE-6 boards of Tokyo University 
Makino et al.EoO^i . Use of GRAPE-6 hardware is essential 
to perform A-body simulations with more than 10 5 stars 
within a reasonable amount of co mputer time. D etails on 
the NBODY4 code can be found inlAarset l JW99|) and ref¬ 
erences therein and IBaumgardt fc MakincT 2003il . 

We first checked the simplest situation, that of a single¬ 
mass cluster. In Fig. 21 we compare the evolution of the La- 
grange radii of a Plumm er model as obtained with NBODY4 
feaumgardt et ai]|2003 1 and ME(SSY)**2. The agreement 
between the two methods is excellent if we set the coefficient 
in the Coulomb logarithm to 7 C = 0.09 when converting the 
A-body time units, natural to NBODY4, to FP units. Given 
the run-to-run variations of MC results (for different ran- 
do m sequences), thi s value is compatible with the one found 
bv lGiersz fc Heg’gie 1 ( 1 1994 !) in comparisons between A-body 
runs with various part ic le numbers ( A p = 250 — 2000) and 
later confirmed bv IBaumgardt! i200l ) (A p = 128 — 16 384), 
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Figure 5. Evolution of Lagrange radii during the core collapse 
of a Plummer model with a mass spectrum. The cluster has a 
Kroupa mass function (see text) extending from 0.1 to 10 Mq. 
We compare the results of ME(SSY)**2 with N p = 10 6 particles 
(solid lines) to those of a direct IV-body run with N p = 131 072 
(dotted lines, in magenta in the on-line colour version). A value of 
7c = 0.015 was used in the Coulomb logarithm when converting 
A-body time units into FP units to get the best overall agreement. 



Figure 6. Mass segregation during the core collapse of a Plum¬ 
mer model with a mass spectrum. The data is for the same MC 
and A-body simulations as in Fig. 0 We plot the stellar mass 
averaged over all particles inside spheres containing the indicated 
fraction of the total mass of the cluster. 


Figure 7. Evolution of Lagrange radii during the core collapse 
of a Plummer model with a broad mass spectrum. The cluster has 
a Kroupa mass function (see text) extending from 0.1 to 100 Mq. 
We compare the results of ME(SSY)**2 with N p = 1.25 x 10 6 
particles (solid lines) to the averaged results of two direct AT-body 
runs with N p = 131072 (dotted lines, in magenta in the on-line 
colour version). Our standard value of 7 C = 0.01 was used in the 
Coulomb logarithm when converting A-body time units into FP 
units. 



Figure 8. Mass segregation during the core collapse of a Plum¬ 
mer model with a broad mass spectrum. The data is for the same 
MC and A-body simulations as in Fig. Q See caption of Fig. IH1 
for explanations about the plotted quantities. 
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7 c = 0.11 ; it is als o compatible with the one chosen by 
iDrukier et alJ dl99fiii to adjust their FP results to A-body 
simulations, 7 c = 0.10. 

After this test-run, we considered a cluster consisting 
of stars with masses from 0.1 to 10 Mg. The masses are dis - 
tributed according to the IMF advocated bvfKrouDal d2QQll) ■ 
For our simulations, the “Kroupa IMF” corresponds to a 
piecewise power law: dA*/dM* oc M* a with a = 1.3 be¬ 
low 0.5 Mg and a = 2.3 for higher masses. This produces a 
mass ratio p ~ 19.3. Fig. 0 depicts the Lagrange radii evo¬ 
lution for ME(SSY)**2 and NBODY4 simulations of this 
model. Converting time units with 7 C = 0.015, we observe 
an excellent agreement between the results of both codes un¬ 
til, at f ~ 0.043 Tfp, shortly before the moment of deepest 
collapse in the A-body run, it starts showing slower con¬ 
traction of the innermost regions. This departure from the 
MC results is probably due to small-number effects, such as 
large-angle scatterings or binary formation, that naturally 
kick in in the shrinking core of the A-body system (com¬ 
puted with N p = 131 072) but are, by nature of the MC 
approach, absent from the ME(SSY)**2 run. Eventually, at 
t ~ 0.047 Tfp, 3- body binaries r everse core collapse in the A- 
body simulation. iHenon £ 1975 1 explained why a smaller -y c 
value is appropriate for systems with a mass spectrum and, 
from multi-mass A-bod y simulations with N p = 250— 1000, 
lOiersz fc Heggiel (Il99fitl indeed found 7 C ~ 0.015 — 0.025. 

In such multi-mass clusters, core collapse is driven by 
mass segregation. Fig. fallows us to witness this process by 
plotting the evolution of the average mass within spheres 
containing various fractions of the total cluster mass. Again 
the MC results match those of the A-body run closely. 

We now consider a mass spectrum of realistic breadth, 
i.e., a Kroupa IMF extending from 0.1 to 100 Mg, yield¬ 
ing p ~ 157. To reduce numerical noise without increasing 
N p to an impractically high value, we realised two A-body 
simulations of this system with A p = 131 072 (starting from 
different realisations of the initial cluster) and averaged their 
results. Comparison of the Lagrange radii and average stel¬ 
lar mass evolutions are shown in Figs Q and |S] This time, 
we stuck to 7c = 0.01, the value we traditionally use to 
convert from FP time to physical time in our MC simula¬ 
tions of multi-mass clusters. This value turns out to yield a 
perfect match between the MC core-collapse time and the 
time of maximum contraction of the inner regions in the 
A-body runs. However, -y c — 0.025 would have allowed a 
better agreement in the early shapes of the Lagrange radii 
curves. At any rate, the agreement between ME(SSY)**2 
and NBODY4 results is not as good as for previous cases. 
In particular, it appears that mass segregation is faster and 
stronger but more progressive in the A-body run. In the 
MC simulation, segregation accelerates at late times and 
the average stellar mass in the innermost regions (inside the 
1 % Lagrange radius) reaches values similar to those found 
in the A-body run near the moment of collapse. A better 
understanding of the cause of the differences between the 
two simulation methods in the regime of broad IMF may 
be reached in future investigations thanks, in particular, to 
A-body simulations with higher A p and, hence, less noise 
and less affected by small-number effects in the central re¬ 
gions. For the moment, we note that the ME(SSY)**2 evo¬ 
lution of this broad-IMF cluster is qualitatively similar to 
that shown by the direct A-body integration and that good 



Figure 9. Evolution of Lagrange radii during the core collapse 
of a Plummer model with a broad mass spectrum. The cluster 
has a Salpeter mass function extending from 0.2 to 120 Mg. We 
compare the results of ME(SSY)**2 with A p = 1.25 X 10® par¬ 
ticles (solid lines) to those obtained with SPEDI (with dotted 
lines, in magenta in the on-line colour version). In order to allow 
a better comparison of the shapes of the curves, we aligned the 
core collapse times by multiplying the time units for the SPEDI 
run by a factor 1.15. 


quantitative agreement is obtained for the aspects most im¬ 
portant to the present investigation of the collisional run¬ 
away, namely the core collapse time and the degree of mass 
segregation reached in the central regions during late col¬ 
lapse. The core collapse time we obtain with ME(SSY)**2 
is tec ~ 7.3 x 10 -3 Tfp = 7.9 x 10 _2 t rh (0) = 0.17f rc (0), 
in agreement with the value of t cc — 0.15t r c(0) found in 
GFR04 for systems with p > 50. 

As for the value of 7 C , we decided to stay on the con¬ 
servative side by keeping it at 0.01. This is probably slightly 
too small, hence predicting a core-collapse evolution too slow 
with respect to other time scales, most importantly that for 
stellar evolution, but the difference in relaxation time com¬ 
pared to, say, 7 c — 0.025 is smaller than 10% for A* ^ 10®. 


4-2.2 Comparison with the gaseous model 

The direct A-body method represents the most accurate 
but most computationally expensive way of simulating the 
secular evolution of a stellar cluster, subject to relaxation. 
Unfortunately, precisely because it treats gravitation in such 
a direct fashion, it offers no or little clean and easy way of 
establishing the global behaviour that a system may exhibit 
in the limit of very large number of stars (A* 10 5 ) from 

the results of simulations made with much fewer particles 
and ridden with noise and small-A effects. 

A nearly opposite approach is taken in SPEDI, a cluster 
evolution code developed by R. Spurzem and coll aborators 
llLouis fc Spurzenilll99ll : iGiersz fc Spurzeinit 19941 : Spurzem 



































Runaway collisions 


15 


1995 ; lAmaro-Seoane et alJl20o3) 7 . Here, like in 


& Takahashi 

the MC scheme, the star system is treated in a explicitly sta¬ 
tistical manner which assumes a very large number of stars. 
By taking velocity moments of the Boltzmann equation up 
to the second order, one obtains a set of partial differen¬ 
tial equations for the evolution of density, average (radial) 
velocity and velocity dispersions (radial and tangential) of 
the stars at each position in the (spherical) cluster. These 
equations are similar to those governing the structure of 
self-gravitating spherical gas cloud. The effects of 2-body 
relaxation are treated through local prescriptions for heat 
conduction among stars of the same mass and energy ex¬ 
change between stars of different masses. This method, al¬ 
though approximative, allows fast simulations and smooth 
results because the cluster is treated as a continuum. In this 
scheme, the mass spectrum has to be discretized into a num¬ 
ber of components, each representing stars of a given mass. 

In Fig. |^| we compare the evolution of Lagrange radii, as 
obtained with ME(SSY)**2 and SPEDI, for a cluster with 


a Salpeter IMF (dN*/d,M* oc M, 


extending from 0.2 


to 120 Mg. The SPEDI simulation is the same as presented 
in Fig. 2 of GFR04; it was carried out using JV CO mp = 50 
mass components. Our tests show that the results depend 
very little on JV comp as soon as more than ~ 12 components 
are used. We have set the gaseous model parameters to their 
standard values: the c oefficient in the thermal conductivity 
relation is A = 0.4977 llGiersz fc Spurzem|l994|) and c ouipar - 
tition time is set by A eq = 1 IffiuimenM^^Tffiahash 1119951. 
We note that SPEDI produces a collapse some 13 % faster 
than ME(SSY)**2, a difference that may be the result of 
an slightly inadequate value of A eq as this parameter has 
only be determined for the case of 2-component models, by 
comparisor^vM^VUao^^m^J^^nmxdaUon^Spurzem & 
Takahashi ll995i^ ^rTh^knmj^w!Hiave^esc£ded the time of 
the gaseous-model simulation to obtain the same t cc as in 
the MC run. The evolution to core collapse is just slightly 
more progressive in the SPEDI run but the agreement is 
otherwise very good, given the important differences in both 
numerical approaches. The core collapse time we obtain is 
tec ^ 5.4x 10 -3 Tfp = 5.8x 10“ * 2 * * * * t rh (0) = 0.12t rc (0), ~ 20% 
faster than in GFR04. This difference reflects not only differ¬ 
ences between the two MC codes but also intrinsic variations 
in tec obtained between ME(SSY)**2 runs computed for the 
same initial conditions but with different random sequences 
(see Fig. 2 of Paper II). 


4.3 Runaway collisions: Comparison with Quinlan 
& Shapiro (1990) 

In iFreitag fc Bend ll2002blf . we have already checked that 

ME(SSY)**2 produces the correct rate of collisions, both as 

a function of the distance to the centre of the cluster and 

of the masses of stars, for multi-mass clusters. We have also 

reproduced the overall results of collisio nal models of M BH- 

hosting g alactic nuclei consid ered by iDuncan fc Shapirol 
dl983h and lMurohv et ah! dl99ll) . In such environments, how- 


7 SPEDI is maintained at the Astronomisches Rechen- 
Institut in Heidelberg, see http://www.ari.uni-heidelberg.de/ 
gaseous-model. 
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Figure 10. Evolution of Lagrange radii in core collapse with 
stellar collisions. The initial conditions are those of model E4A 
of QS90. We present two simulations. In the first one (solid lines, 
QS90-E4A), all collisions are treated as pure mergers with complete 
rejuvenation. In the second case (dashed lines, QS90-E4Asph), we 
use SPH-inspired collision prescriptions and minimal rejuvena¬ 
tion. 


ever, collisions occur at very high relative veloci ties and do 
not lead to runaway gr owth lFYeitag_et_alJ 200jh 

To our knowledge, iQuinlarT^^ShapiTO l jl99(l hereafter 
QS90), using a FP code, were the first to study in a sys¬ 
tematic way the process of collisional runaway accounting 
self-consistently for cluster dynamics and secular evolution 
and how collisions themselves affect it. 

We have run simulations for all 6 “E” models of QS90, 
assuming, in most cases, sticky-sphere collisions and com¬ 
plete rejuvenation for each merger as these authors did. Ini¬ 
tial conditions are high-density Plummer clusters with all 
stars of solar type, with various total masses and sizes. Ta¬ 
ble |21 lists the characteristics of these models. Stellar evolu¬ 
tion is included. QS90 implemented it in a simplified way, 
similar to ours. However, while we can assign an individ¬ 
ual age to each MC particle (which represents 4.2 to 36 
stars, depending on the model), in the FP scheme of QS90 
only a small number of discrete mass component are present 
(for stars of mass 2* Mq with i = 1, 2,... 8), each of which 
represents some homogeneous population. Stellar evolution 
can only be treated in statistical way, with each component 
having an average age. Similarly, collisions can only be ac¬ 
counted for statistically by integrating merger rates within 
and between mass components and distributing the number 
of merger products over components in a way that conserve 
total mass but not necessarily individual stellar masses be¬ 
cause not all possible merger masses are represented. Unlike 
QS90, we have allowed for collisions between MS stars and 
compact remnants by assuming complete disruption of the 
MS star but no change to the mass or velocity of the compact 
object. Certainly an oversimplification, this prescription ob¬ 
viously corresponds to the most unfavourable possible one, 
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Table 2. Cluster simulations for comparison with lQuinlan &; ShaoirJ il99Glh 


Name 

Name in QS90 

TV* 

TV P 

T?nb 

(pc) 

tec 

(Tfp) 

tcc 

(Myr) 

QS90-E4A 

E4A 

1.8 x 10 7 

5 x 10 s 

0.41 

0.134 

141 

QS90-E4Ahr 

E4A 

1.8 x 10 7 

2.1 x 10 6 

0.41 

0.129 

136 

QS90-E4Asph( a ) 

E4A 

1.8 x 10 7 

5 x 10 5 

0.41 

0.275 

289 

QS90-E4B 

E4B 

3.1 x 10 7 

5 x 10 s 

0.71 

0.182 

553 

QS90-E4Bb (b ) 

E4B 

3.1 x 10 7 

5 x 10 s 

0.71 

0.172 

522 

qS90-E4Ac( c > 

E4B 

3.1 x 10 7 

5 x 10 5 

0.71 

0.189 

574 

qS90-E2A 

E2A 

6.2 x 10 6 

5 x 10 5 

0.565 

0.371 

397 

qS90-E2B 

E2B 

1.1 x 10 7 

5 x 10 s 

1.00 

(d) 


qS90-ElA 

E1A 

2.1 x 10 6 

5 x 10 5 

0.765 

0.758 

803 

qS90-ElB 

E1B 

3.6 x 10 6 

5 x 10 s 

1.31 

1.23 

3680 

qS90-E4Ab( e ) 

E1B 

3.6 x 10 6 

5 x 10 5 

1.31 

1.24 

3710 


Collision prescriptions based on SPH data. Minimal rejuvenation. 
Other random sequence than for QS90-E4B. 

( c ) Time step four times larger than for QS90-E4B. 

( d ) Core collapse stopped by stellar evolution. 

Other random sequence than for QS90-E1B. 



Figure 11. Number of stars in various mass bins during cluster 
evolution with stellar collisions. Same simulation as in Fig. cm 
(QS90-E4A). Dotted lines are the results of QS90 (their Fig. 2a) for 
mass-component of mass 1, 2, 4, 8 and ^ 16 M®. The curves for 1 
and 2Mq are nearly indistinguishable from ours (for 1 — 1.4 and 
2 — 3Mq bins). We obtain a lower number of stars more massive 
than 3 Mq than QS90. At the end of our simulation, ~ 200 white 
dwarfs and a similar number of neutron stars have formed. The 
(negligible) contribution of neutron stars is included in the first 
mass bin; white dwarfs are not plotted. We used 5 x 10 5 particles 
for our simulation so each of them represents 36 stars. 


as far as collisional runaway is concerned. To have the same 
ratios between the time scale for relaxation and that for 
other processes (collisions, stellar evolution) as QS90, we 
adopt their 7 C = 0.4 value for this set of simulations. 

In contrast to QS90 we find core collapse and collisional 



Figure 12. Mass density of stars of various masses in the central 
region of a cluster during its evolution with stellar collisions. Same 
simulation as in Figs [IT] and rmroS9Q-E4Ai. For the 1 - 1.4Mq, 
2 — 3Mq, 4 — 7Mq, 8 — 15 Mq and 16 Mq mass bin, we plot 
the average density within a sphere containing, respectively, a 
fraction 0.001, 0.01, 0.01, 0.1 and 0.1 of the total mass in the 
corresponding range of stellar masses. As in Fig. El thin dotted 
lines indicate results of QS90 (their Fig. 2b). Note that the curve 
from QS90 closest to our line for the 4 — 7Mq bin is for their 
8Mq mass component. 

runaways in all cases but E2B whose core contraction is 
stopped by the mass loss due to stellar evolution of collision 
products. 

Agreement for model E4A is very satisfying: in our 
model QS90-E4A, we find core collapse and runaway to 
happen at t cc — 140Myr ~ 1.44f r h(0), to compare with 
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Figure 13. Number of stars in various mass bins during cluster 
evolution with stellar collisions. Our simulation for model E4B 
(QS90-E4B) is compared to data from QS90 (their Fig. 4a). See 
caption of Fig. [lljfor explanations. 

QS90’s value of ~ 130 Myr. Without collisions, a single¬ 
mass Plummer model would experience core collapse at 
t cc — 17— 18t r h(0) so mergers clearly play an important role 
in the cluster evolution from the beginning, by dissipating 
energy and creating a mass spectrum, hence allowing mass 
segregation. For this reason, one cannot predict the occur¬ 
rence of runaway in those clusters just from the relaxational 
value of the core-collapse time, £ cc | rlx . Also, even though all 
stars are initially of one solar mass, it is not sufficient to have 
£cc| rlx < 10 Gyr for collisional runaway to kick in. Paradox¬ 
ically, if merger products are formed before deep collapse 
they may be able to stop collapse as they evolve off the MS 
much earlier, as happens in our case E2B (QS90-E2B). 

Fig. ED] shows the evolution of the Lagrange radii our 
simulation of case E4A. We have also redone the simulation 
with the more realistic prescriptions for collision outcomes 
based on SPH simulations (QS90-E4Asph). As can be seen in 
Fig. ESI when collisional mass-loss and, more importantly, 
fly-by, non-merging collisions are accounted for, the cluster 
evolution is significantly slower. This is because, in clusters 
with such a high velocity dispersion (ct„, c ( 0) ~ 231 kms -1 ), 
only a minority of collisions result in mergers, amounting to 
a reduction in the effective collision cross sections. In Fig. EH 
we have plotted the evolution of the number of stars of vari¬ 
ous masses during the evolution of QS90-E4A. The evolution 
of the number of stars of mass 2 to 3 M@ in our run is nearly 
identical to that of the 2 Mq component of QS90 but we ob¬ 
tain significantly fewer stars more massive than 3Mq. This 
demonstrates that the rate of collision between 1 Mq (and 
hence, of formation of first-generation mergers) is the same 
in both simulations. Because QS90’s FP method imposes a 
completely different (and much less physical) way of dealing 
with mergers of mergers, it is not surprising no close agree¬ 
ment is reached on that matter. The M-R relation cannot be 


blamed, however, as they assume R,/ Rq = (M*/Mq) 0 ' 55 , 
which is close to our relation. The drop in the number of 
stars more massive than 16 Mq at the end of our simulation 
corresponds to the run-away growth of a VMS by merging 
between massive stars. When we stopped our run, the VMS 
had reached 651 Mq. We refer to Paper II for a discussion 
of physical processes that may terminate the VMS growth, 
most of which are lacking from ME(SSY)**2. How mergers 
create a population of massive stars which come to domi¬ 
nate the central region of the cluster is further indicated by 
Fig- EH Here, we follow the contribution to the central den¬ 
sity of stars in various mass ranges. Unlike QS90, we cannot 
determine the density at R = 0 but have to sum the mass 
in some small spherical volume to estimate it. Our values 
are therefore lower limits. Still, the agreement for the lowest 
two mass bins (1 — 1.4 Mq and 2 — 3 Mq) is very satisfying, 
but again, we get considerably fewer higher-mass stars. 

Although 3 times less massive, model E4B has the same 
velocity dispersion as E4A and, hence a very similar initial 
value of frix/icoii (up to small differences in Coulomb log¬ 
arithms). Therefore, one could have expected collisions to 
play the same role and to get the same value of £ C c/£rh(0). 
However, stellar evolution introduce still another time scale 
in the problem and because f r h(0) is longer (in years) for 
model E4A, more collision products have time to evolve 
off the MS, which delays core collapse. At t = 0.13 Tfp, 
the fractional number of white dwarfs in our E4A run is 
~ 10 -5 while it reaches ~ 10 -3 in run E4B. In QS90’s run 
for this case, stellar evolution was able to stop core collapse 
around t as 700 Myr, before any star more massive than 
32 Mq formed. In contrast, our models produced a runaway 
merger after 520 to 570 Myr already. The evolution of num¬ 
ber of stars of different masses is plotted in Fig. EH and 
compared with QS90’s results (for the contraction phase). 
Again, the agreement with QS90 is excellent for the 2 —3 Mq 
mass bin but poorer for others. This time, our run produces 
more massive stars (after t as 300 Myr), which eventually 
leads to runaway. Clearly, because the evolution is driven 
by a positive feedback loop between mass segregation and 
collisions and because only stellar evolution of merger prod¬ 
ucts themselves can revert it, differences in the treatment of 
collisions may have very important effects. 

For their other “E” models, QS90 only indicate the end 
state in a very succinct way, giving no details about the 
evolution. For model E2A, we find runaway at t ~ 400 Myr, 
to be compared with QS90’s value of t ~ 490 Myr, a sat¬ 
isfying agreement, given how different the approaches are. 
Model E2B is prevented from entering the run-away phase 
by stellar evolution in both QS90’s and our simulation. In 
our simulation (QS90-E2B), a few particles grow to 8Mq 
and one to 11 Mq during the phase of central concentration. 
QS90 report that the most massive stars formed have 8 Mq 
in this case. In QS90, the core collapse of model E1A is re¬ 
versed by the formation of, and energy input from 3-body 
binaries, a process we do no take into consideration (see dis¬ 
cussion in Appendix 01, but the fact that they form stars 
as massive as 64 Mq indicates that they nearly reach the 
conditions for runaway. For this cluster we obtain runaway 
at £ ~ 800 Myr. Finally, E1B is another case for which QS90 
find the core collapse is interrupted by stellar evolution but 
our simulation produces runaway at t ~ 0.37 Gyr. In one of 
our simulations for this model, the runaway star is destroyed 
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in a collision with a compact object, just after it has reached 
a mass of 157 Mg. Another runaway sequence starts shortly 
afterwards. 

Since QS90, runaway collisions in clusters have been 
studied by usf^rf^if^direc^JVNmc^^mRioc^^Portegies 
Zwart et al. TM99an Portegie^Zw£^^^^lc^nMi|2002^ ^orte- 
gies Zwart et al. 12004) . Unfortunately, these simulations 
were either limited to A p < 10 5 , in which case binaries play 
a crucial role or, for the few cases reaching N p ~ 6 x 10 5 
and experiencing runaway, done for clusters with high initial 
concentration (King parameter Wo sS 9) and, hence proba¬ 
bly dominated by small-number effects in the innermost re¬ 
gions. In Paper II, we report briefly on the simulations we 
have don e with initial conditions sim ilar to those of some 
models of lPortegies Zwart et a,l ] (12004 . 

5 SUMMARY 

The results presented here are mostly tests to ensure that 
the Monte Carlo stella r dynam ics code we use, ME(SSY)**2 
iFreitag fc BendfeOQll 12002b 1. correctly treats the key pro¬ 
cesses at play in the runaway route, which are: 

(i) Mass segregation-induced core collapse, driven by 2- 
body relaxation in cluster with a broad, realistic mass func¬ 
tion (Salpeter of Kroupa, with M»,min = 0.1 — 0.2 Mg and 
M*, max ~ 100 Mg) 

(ii) Effects of collisions in the evolution of the cluster and 
occurrence of collisional runaway. 

T his paper is a companion to Paper II llFreitae et alJ 
l2005h in which we present the results of our large set of 
simulations to determine the conditions for, and character¬ 
istics of collisional runaway in young stellar clusters. In the 
scenario we investigate, the collisional phase is brought up 
by the concentration, through relaxation, of massive stars 
in the centre of a cluster. A star much more massive than 
100 Mg, formed in a quick sequence of collisions may not 
only be a progenitor for an IMBH but is in itself an exotic 
object of considerable interest. 

To address point |(i)| we performed simulations of clus¬ 
ters with ME(SSY)**2, NBODY4 and the gaseous-model 
code SPEDI. The only physical process included in the MC 
and gaseous-model runs was 2-body relaxation, treated in 
the standard Fokker-Planck approximation, i.e., as sum of 
uncorrelated small-angle two-body scatterings whose rate is 
determined by local conditions at each point in the cluster. 
The A-body code treats Newtonian gravitation in an es¬ 
sential exact way, without assumptions about the nature of 
relaxation. ME(SSY)**2 produces results in close agreement 
with their A-body counterparts for fi = M* jmax /(M*) ^ 20 
at least. Realistic IMF correspond to /r > 100, however. 
There are clear discrepancies between the MC and A-body 
simulations in this regime. Most noticeably A-body results 
show an initially stronger but more progressive concentra¬ 
tion of massive stars in the central regions. Nevertheless, the 
most important characteristics of the core collapse as a path 
to the collisional runaway stage, namely the time for it to 
happen, t cc , and the magnitude of central mass segregation 
during deep collapse are similar in both type of simulations. 
Furthermore, SPEDI yields results very close to those of 
ME(SSY)**2, except for a value of t c c 13% shorter. 


Turning now to point |(ii)| we note that only very few 
numerical simulations of the evolution of clusters subject to 
relaxation and collisions, up to and including the runaway 
stage have been publi shed. Putti n g aside the recent A -body 
runaway simulations JPo^gie^^w^^^^aL|^999alPorte- 
gies Zwart & McMilla,n* l2001lTPo^tde^^wmt^^djl2f)f)4 . 
not suitable for clear-cut comparison with ME(SSY)**2 re¬ 
sults as small-A effects may play a stro ng role in them, we 
are lef t with the older FP simulations of lQninla.n fc Shanirol 
( 199(1 QS90). These models lack realism because they start 
with a single-mass population of 1 Mg stars. Collisions are 
assumed to result in perfect mergers, with no mass loss. In 
most cases considered, the cluster is dense enough to pro¬ 
mote stellar mergers at early times. The more massive stars 
thus formed accelerate collapse by mass segregation, hence 
further increasing the merger rate. These same merger prod¬ 
ucts can also terminate the process of central density build¬ 
up through their mass loss when they evolve off the MS, well 
before the original 1 Mg would have done so. 

We have simulated the 6 models of QS90 for which 
these authors have assumed that the gas lost through stellar 
evolution escapes the cluster (rather than staying in it and 
forming new stars). Although collisional runaway happens 
more often in our simulations than in QS90 runs, we ob¬ 
tain good general agreement with the relatively scarce data 
published by QS90. For the two situations in which they 
report core collapse uninterrupted by stellar evolution or 3- 
body binaries, we obtain core collapse times 8 % longer and 
18 % shorter than theirs, which we find satisfying consider¬ 
ing the widely different numerical methods and, especially, 
treatments of collisions, of critical importance here. The pro¬ 
duction rate of 2 — 3 Mg objects in our runs is nearly exactly 
the same as for QS90’s 2 Mg mass bin but noticeable dif¬ 
ferences exist for higher-mass objects, which, resulting from 
longer merger sequences, are more affected by differences in 
treatments of collisions. For a large part, the discrepancies 
between ours and QS90’s certainly originates in this. Our 
treatment of collisions being much more direct and accurate 
than the one allowed by the FP code of QS90, the differences 
found in these comparisons do not cast doubt on the ability 
of ME(SSY)**2 to deal with collisional cluster evolution. 

We have also re-simulated one of the QS90 models with 
the highest velocity dispersion using our more realistic treat¬ 
ment of^ollision^^ase^^i^^PE^imxdaBon^^Freitag & 
Benz l200fih and found a significantly longer core-collapse 
time due to the collisions being less effective at producing 
higher-mass stars. This stresses the importance of using col¬ 
lision prescriptions which account for fly-bys and mass loss, 
in the high-velocity regime. This question is investigated in 
more depth in Paper II. 


APPENDIX A: NEGLECT OF BINARIES 

In principle, “hard” binary stars, either primordial or formed 
during core collapse through 3-body processes, may play 
an important role in the cluster evolution. During gravi¬ 
tational encounters with single stars, a hard binary is likely 
to shrink, thus allowing the single star (possibly a former 
member of the binary if an exchange has occurred) to emerge 
with increased kinetic energy. Through this dynamical heat¬ 
ing, hard binaries may suspend or even reverse core collapse 
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Figure Al. Core collapse of a Wo = 3 King cluster with Salpeter 
mass function, computed with SPEDI. The mass spectrum is dis¬ 
cretized into 15 components with the indicated individual stellar 
masses (m) and mass fraction (/ m ). On the top panel, we repre¬ 
sent the contribution of each mass component to the central mass 
density. The dotted line is the total density. The bottom panel 
shows the velocity dispersion of each component. The dotted line 
is the mass-averaged dispersion. 


llHeggie fc HuJ2003l and references therein). Obviously, this 
might prevent the core from ever entering the high-density 
collisional phase needed in the runaway scenario. However, 
when the finite size of stars is taken into account in the 
numerical study of single-binary and binary-binary interac¬ 
tions, it appears that th e col li sion of at leas t two of the 
stars is a likely event (see lFregeau et aJll2004l for the most 
recent cross-section computations and references about this 
questionb^Hiougl^eing^oin^nas^^iamic^Giersz & 
Spurzem (2003) studied the statistics of binary-single and 
binary-binary encounters occurring in their simulations of 
clusters containing primordial binaries and found that of or¬ 
der 50 % of these interactions should indeed lead to stellar 
mergers, for typical globular cluster parameters. It is there¬ 
fore possible that binaries actually foster collisions rather 
than preventing them. This is exactly what was found in N- 
body simulations of relatively small clusters (most of them 
with IV* ^ 131 072, one run with 77* ~ 585 000) by Porte- 
gies Zwart and collaborators d Portegies_Zwart_et ~aflll999bt 

Portegies Zwart fc McMillanll2002nPortedei^Zwari^^dJ 

^O^b Furthermore. dissipative processes, such as collisions 
or tidal interactions, occurring during binary interactions 
may significantly reduce the heating effect of hard bina- 
ries^omDare^^iH^h^£oin^nus^£aroximatioi^McMil- 

Ian TM8ftlGoo^ma^j^^nimds^T^9dnChenio^fc Huang! 
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Figure A2. Same gaseous model simulation as on Fig. cm Here 
we plot, for each mass component, the central collision rate per 
star (top panel) and integrated number of collision or collision 
probability (bottom panel) for a star near the centre. To avoid 
double counting, the rate for component i includes collisions with 
all components j i. The steep (black) dash-dotted line is an 
estimate of the formation rate of 3-body binaries (per star) ob¬ 
tained by using average stellar density and mass in equation cm . 
Note that this simulation do not include the effects of collisions 
or binary formation but only 2-body relaxation. Collision and bi¬ 
nary formation rates have been estimated a posteriori assuming 
the cluster is made of 10® stars and its size is 7 ?nb = 0.2 pc. 

Il996tl . To our knowledge, however, the impact of this on core 
collapse has not yet been studied explicitly. 

Even if binaries probably do not prevent core collapse 
and collisions, it does not follow that they will not im¬ 
pede or modify the runaway process. Indeed, the cross sec¬ 
tion for collision of a single star with a binary is of order 
nGaM(Vf^)~ 2 oc M where a is the binary semi-major axis 
and M the mass of the three stars, while it is approximately 
ttGRM(V ^)~ 2 oc M a with a ~ 1.5 for the collision with 
a single more massive star of mass M (> 50 Mq) and ra¬ 
dius R. According to mathematical modelling through the 
coagulation equation, runaway is expected only if the cross 
section scalefGike_^^^iH^^^ 1 iLeJl99.1l200nt Malyshkin 
& Goodman 1200 ill .This condition is not obeyed in the case 
of binaries competing with each other for interaction and 
merger with single stars. However, the coagulation equation 
is clearly at best a crude idealisation for the complex stellar 
dynamical situation of interest here. 

To simplify the problem, we have assumed in this work 
that no primordial binaries were present. We now examine 
if we are then justified to neglect binary process altogether 
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consistently, in spite of the possibility of forming binaries 
dynamically. We do not consider binary formed by tidal in¬ 
teractions because their formation requires close passage at 
a distance of a few stellar radii at most. Their formation rate 
cannot then be significantly higher than collision rate and 
their post-capture evolution, though still a controversial is¬ 
sue, is likely to lead to merger anyway. We instead consider 
the question of binaries formed through (point-mass) 3-body 
interactions. 

The c reatio n rate of such binaries is ( Binney fe 
Trem aine TTosT ^q - f8^ir see also appendb^f l^mova et ahl 


fl3b — C3bri' 


G 5 Mi> 


(Al) 


where n is t h^^^la^^^^r^nsi^^ i^^lg^ 
We compare this with the collision rate, 


: 0.75 (Good- 


^-coll 


One gets 


nt 00 u ~ Rt ^1 + 
n 2 GM,R * 


2 n 2 


2 GM, 

R*(V%) 2 




87T- 


flcoll _ 87T ^ 500 /<T„\ 

~ C 3b G 4 A-/*n ~ Fm [v*J 


(A2) 


(A3) 


with Vq = ( 2 GM 0 /R 0) 1 / 2 = 618 km s x . For typical 
100 M 0 stars (dominating the central regions), f?* ~ 14R 0 , 
V, ~ 1670 km s" 1 and 


ficoll 

n 3 b 


100 M© 


10 6 pc - 


20 km s 


(A4) 


Thus, one cannot clearly exclude that the formation of 3- 
body binaries will not compete with direct collision, at least 
for systems with a relatively low velocity dispersion. 

Estimating a posteriori the central values of n^h from 
MC runs is very difficult because of the steep dependences on 
n and a v , which make the estimate extremely noisy. Thus, 
to see how how h co \\ and 723 b evolve during core collapse, 
we resort to the SPEDI gas-model code presented in Sec¬ 
tion rrm In Figure. EH we follow the evolution of cen¬ 
tral densities and velocity dispersions in a 15-component 
SPEDI core-collapse simulation. Assuming N, = 10 6 and 
_R nb = 0.2 pc, we can now compute what the central colli¬ 
sion and 3-body formation rates would have been during the 
core collapse. Because the evolution speeds up near the mo¬ 
ment of core collapse, we use the central potential instead of 
time as independent variable in Fig. EH where we plot the 
instantaneous and time-integrated rates for all mass com¬ 
ponents. For this size and star number, the cluster should 
become collisional before the first binary forms. How this 
depends on the cluster parameters is expressed by either of 
the scalings 


oc nAr ^ 1 oc (In AT rh )“ 2/3 TV, 10/3 (A5) 

723 b 

where R c 1 is some characteristic cluster size (e.g., the half¬ 
mass radius). The strong dependence on A r , suggests that 
3-body binaries, which appear to dominate the collision pro¬ 
cess in TV-body simulations may be of little importance in 
larger systems with TV* > 10®, such as very massive young 


clusters or proto-galacti c nuclei. This analysis is^n agree¬ 
ment with the results of lPortegies Zwart et alJ ( 2004 1 who 
find the dynamically-formed binaries to play a lesser role 
in comparison with^rnffiei^nTOi^^imulaUnns^Portegies 
Zwart & McMillan TlOffi?^ ^ 
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